SpeakEasy code from Gaiteri et al, 2015.
library(tidyverse)
library(janitor)
library(readxl)
library(WGCNA)
net_dir = "/pastel/projects/speakeasy_dlpfc/SpeakEasy_net_MF/"
expression_dir = "/pastel/projects/speakeasy_dlpfc/"
dir_mod_02 = "/pastel/projects/speakeasy_dlpfc/SpeakEasy_net_MF/"In order to summarize the SE network information we will perform the module eigengene approach. Module eigengene is defined as the first principal component of the expression matrix of the corresponding module.
## Read the metabolites matrix and SE results
load(paste0(expression_dir, "exprdata_byregion.Rdata"))
# Expression data for a single set in the form of a data frame where rows are samples and columns are genes (probes):
expr_matx_t = as.data.frame( t(exprData_MF)) # Residuals of the expression
k_dataset_02 = read.table(paste0(dir_mod_02, "geneBycluster.txt"), header = T, stringsAsFactors = F) # clusters from SpeakEasy
k_dataset_ord = k_dataset_02[match(colnames(expr_matx_t), k_dataset_02$ensembl), ] #order the matrices
# all(colnames(expr_matx_t) == k_dataset_ord$ensembl) # must be true Showing the module eigengenes in a dataframe, with each column corresponding to one eigengene.
colors = k_dataset_ord$cluster_lv3
lv3_moduleEigengenes = moduleEigengenes(expr_matx_t, colors, verbose = 0)
# save results
save(lv3_moduleEigengenes, file = paste0(net_dir, "lv3_moduleEigengenes.Rdata"))
write.table(lv3_moduleEigengenes$eigengenes, file = paste0(net_dir, "lv3_moduleEigengenes.txt"), sep = "\t", quote = F, row.names = T)
createDT(lv3_moduleEigengenes$eigengenes)A dataframe in which each column corresponds to a module, with the component varExplained[PC, module] giving the variance of module explained by the principal component no. PC. The calculation is exact irrespective of the number of computed principal components. At most 10 variance explained values are recorded in this dataframe.
How much of the 1st PC explains expression of the module? For example, the Module 01 eigengene explains 15% of the expression of the genes in the module.
Eigengene of a “Module” with one gene is equal 1.
write.table(lv3_moduleEigengenes$varExplained, file = paste0(net_dir, "lv3_variancer_exp.txt"), sep = "\t", quote = F, row.names = T)
colnames(lv3_moduleEigengenes$varExplained) = gsub("X", "M", colnames(lv3_moduleEigengenes$varExplained))
createDT(lv3_moduleEigengenes$varExplained)A dataframe containing average normalized expression in each module.
# The columns are named by the corresponding color with an "AE" prepended, e.g., AEturquoise etc.
colnames(lv3_moduleEigengenes$averageExpr) = gsub("AE", "AE_M", colnames(lv3_moduleEigengenes$averageExpr))
write.table(lv3_moduleEigengenes$averageExpr, file = paste0(net_dir, "lv3_averageExpr.txt"), sep = "\t", quote = F, row.names = T)
createDT(lv3_moduleEigengenes$averageExpr)sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 8
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.12.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] WGCNA_1.70-3 fastcluster_1.2.3 dynamicTreeCut_1.63-1
## [4] readxl_1.3.1 janitor_2.1.0 forcats_0.5.1
## [7] stringr_1.4.0 dplyr_1.0.8 purrr_0.3.4
## [10] readr_2.1.2 tidyr_1.2.0 tibble_3.1.6
## [13] ggplot2_3.3.5 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-3 ellipsis_0.3.2 snakecase_0.11.0
## [4] htmlTable_2.4.0 XVector_0.34.0 base64enc_0.1-3
## [7] fs_1.5.2 rstudioapi_0.13 DT_0.20
## [10] bit64_4.0.5 AnnotationDbi_1.56.2 fansi_1.0.2
## [13] lubridate_1.8.0 xml2_1.3.3 codetools_0.2-18
## [16] splines_4.1.2 doParallel_1.0.17 cachem_1.0.6
## [19] impute_1.68.0 knitr_1.37 Formula_1.2-4
## [22] jsonlite_1.7.3 broom_0.7.12 cluster_2.1.2
## [25] GO.db_3.14.0 dbplyr_2.1.1 png_0.1-7
## [28] compiler_4.1.2 httr_1.4.2 backports_1.4.1
## [31] assertthat_0.2.1 Matrix_1.3-4 fastmap_1.1.0
## [34] cli_3.2.0 htmltools_0.5.2 tools_4.1.2
## [37] gtable_0.3.0 glue_1.6.1 GenomeInfoDbData_1.2.7
## [40] Rcpp_1.0.8 Biobase_2.54.0 cellranger_1.1.0
## [43] jquerylib_0.1.4 vctrs_0.3.8 Biostrings_2.62.0
## [46] preprocessCore_1.56.0 crosstalk_1.2.0 iterators_1.0.14
## [49] xfun_0.29 rvest_1.0.2 lifecycle_1.0.1
## [52] zlibbioc_1.40.0 scales_1.1.1 hms_1.1.1
## [55] parallel_4.1.2 RColorBrewer_1.1-2 yaml_2.3.5
## [58] memoise_2.0.1 gridExtra_2.3 sass_0.4.0
## [61] rpart_4.1-15 latticeExtra_0.6-29 stringi_1.7.6
## [64] RSQLite_2.2.10 S4Vectors_0.32.3 foreach_1.5.2
## [67] checkmate_2.0.0 BiocGenerics_0.40.0 GenomeInfoDb_1.30.1
## [70] rlang_1.0.1 pkgconfig_2.0.3 bitops_1.0-7
## [73] matrixStats_0.61.0 evaluate_0.15 lattice_0.20-45
## [76] htmlwidgets_1.5.4 bit_4.0.4 tidyselect_1.1.2
## [79] magrittr_2.0.2 R6_2.5.1 IRanges_2.28.0
## [82] generics_0.1.2 Hmisc_4.6-0 DBI_1.1.2
## [85] pillar_1.7.0 haven_2.4.3 foreign_0.8-81
## [88] withr_2.4.3 survival_3.2-13 KEGGREST_1.34.0
## [91] RCurl_1.98-1.6 nnet_7.3-16 modelr_0.1.8
## [94] crayon_1.5.0 utf8_1.2.2 tzdb_0.2.0
## [97] rmarkdown_2.11 jpeg_0.1-9 grid_4.1.2
## [100] data.table_1.14.2 blob_1.2.2 reprex_2.0.1
## [103] digest_0.6.29 stats4_4.1.2 munsell_0.5.0
## [106] bslib_0.3.1